#include<iostream>
using namespace std;
#include<cmath>
#include<string>
#include<fstream>

double RK4_step(double xi, double h, double (*func)(double x, double y), double yi){
	double K1, K2, K3, K4;
	K1 = func(xi, yi);
	K2 = func(xi+h/2, yi+h/2*K1);
	K3 = func(xi+h/2, yi+h/2*K2);
	K4 = func(xi+h, yi+h*K3);
	return yi + h/6*(K1+2*K2+2*K3+K4);
}

/*
 * harmonic oscillator: 
 * 	x'' = -x,
 * we can define v = x', v' = -x,
 */

double 

double f(double x, double y){
	return 10 - y*y;
}

int main(){

	int i,n=10000;
	double a=0, b=2, h;
	double *y = new double [n+1];
	RK2(a, b, n, f, 0, y);

	h = (b-a)/n;
	ofstream fp("Raindrop.txt");
	for(i=0;i<=n;i++)
		fp<< a+i*h << "\t" << y[i]<<endl;
	fp.close();

	delete [] y;
	return 0;
}
